N94- 35606 


Modal Decomposition of 
Hamiltonian Variational Equations 

William E. Wiesel 

Department of Aeronautics and Astronautics 
Air Force Institute of Technology 
2950 P Street 

Wright - Patterson AFB > Ohio 4 S 4 SS 


Abstract 

Over any finite arc of trajectory, the variational equations of a Hamiltonian system can be separated 
into “normal” modes. This transformation is canonical, and the Lyapunov exponents over the trajectory 
arc occur as positive / negative pairs for conjugate modes, while the modal vectors remain unit vectors. 
This decomposition effectively solves the variational equations for any canonical, linear time-dependent 
system. As an example, we study the Voyager I trajectory. In an interplanetary flyby, some of the 
modal variables increase by very large multiplicative factors, but this means that their conjugate modal 
variables decrease by those same very large multiplicative factors. Maneuver strategies for this case are 
explored, and the minimum Av maneuver is found. 


1 Introduction 


A Hamiltonian dynamical system can be written as a vector set of differential equations 


X = 


z 


9H 

ax’ 


where X T = (g,-, p,) is termed the state vector, and the matrix Z is 


( 1 ) 


Z = {-'«}' W 

Introduce the small displacement x(f) = X(t) - X 0 (<) from a known trajectory X 0 (<). Then, to first order 
in small quantities, the displacement vector obeys the variational equations 


x = A(t)x = 


,d 2 H 


8X 2 


X 0 


( 3 ) 


As a set of linear equations, the variational equations are formally solved by the fundamental matrix $(t, to) , 
which satisfies 

$ = A(t')*, *(f 0 ,<o) = /. (4) 

Then, the general solution to (3) can be written as x(t) = Q(t, to)x(to). 


2 The Modal Transformation 

In this section we will review the recent discovery of the modal transformation for general time dependent 
linear systems, Wiesel [6], and we will establish this transformation as a canonical change of coordinates. 


m*:**/^ rage blank not filwcd 
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The stability of a general trajectory of a linear system is determined by the Lyapunov exponents. These 
" etl * v,lu “ . _ i , , ,5) 

k(i 0 )i ' u 

extremalized over all initial displacements x<(<o)- Usually (5) includes a limit as tj — ► oo, but not here. 
The restriction to finite time intervals is an absolute necessity, since prediction of chaotic systems is only 
possible for a finite time interval. 

We wish to find the vectors x<(to), which extremalize the growth of the norm of displacement vectors, 
| x (t,)| with respect to the initial displacement x(to). This is a constrained maximization, since in a linear 
system we may specify |x(i 0 )| = 1 from the outset. Using a Lagrange multiplier n, we have the optimization 

problem . 

J — |x(ty ) | 3 — /i (|x(to )| 3 — l) • («) 

Now, since x(f/) = $(f/»fo) x (*o)> scalar function (6) becomes 




Partial derivatives can now be calculated as if all components of the initial conditions z 0 < were independent, 
yielding 

- ^ = 0 = V x °> ~ t 1 ** 0 * ’ w 

2 dx 0 k Y j 

where *=1,2, ....N. But this is just the component form of 

{$ T $ - ml) e<(<o) = 0. ( 9 ) 

That is, the e<(to) are the real, orthogonal eigenvectors of the real symmetric matrix $ T $, or the right 
singular vectors of <&. Comparison to (5) shows that the Lyapunov exponents over the time interval (to, tj) 

are found from , . 

Hi = exp {2A,(t/ - to)} ■ ( 10 ) 

This has been recognized by Goldhirsch, Sulem and Orszag [1]. We will refer to our A,- as regional Lyapunov 

exponents, since they pertain to the finite time interval (to, t/). 

A matrix $ is symplectic if it obeys = Z, or equivalently = Z. It is well known that the 

fundamental matrix $ is symplectic for a Hamiltonian dynamical system, see, e.g., Wiesel and Pon en [ j. 
But then examining $ T $, we find 

= ^ T Zi = Z, ( 11 ) 

so that is itself symplectic. The eigenvalues of a symplectic matrix occur as inverse pairs, m, 1 /#*<> 90 
by (10) the regional Lyapunov exponents occur as positive / negative pairs. Since the Lyapunov exponents 
are also real, at most half of the modes are unstable, while the other half are stable. The proof of Liouville s 

theorem follows from this as a very simple consequence. 

Over a finite arc of the trajectory, the regional Lyapunov exponents may be used to factor the dynamics 
into separate modes. The initial conditions e,(t 0 ) introduce N special solutions to the variational equations, 
x,(t) = $(< to)e,(<o), on which the average exponential rate of expansion or contraction is an extremum. 
But local variations in these rates can be quite large, Haubs and Haken [2], Nese [3], Sepulveda, Badii, and 
Poliak [4]. We wish to use these N special solutions to the variational equations as basis vectors for the 
entire solution set, and it would be very inconvenient for them to be anything other than unit vectors. Their 
instantaneous rates of change of magnitude are given by 


*<(0 = 


Xj • Ax,- 
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Since the regional Lyapunov exponents are the average of these instantaneous rates on these N extremal 
solutions, we have 

«> 

Then, define N new functions ©*(*) as the solutions to 

©•(*) = Aei - (14) 

with initial conditions ei(fo) on the interval (to,t/)« They are, by (14), trivially unit vectors on the entire 
interval to < t < tj, and orthonormal at t = to. That they must also be orthogonal at t = tj can be seen 
by realizing that e,(t/) must also be the extremal initial conditions for exponential growth of trajectories 
running backwards in time. So they are the eigenvectors of the symmetric matrix (#“ 1 ) T ( < J“ 1 ), and are 
orthogonal. But at other times in the interval ( to, if ) the e,*(i) vectors may not be orthogonal. We note 
that since the new vectors remain unit vectors, 

<Ti(f) = ©i • Ae,* (15) 

is an alternate form of (12). The e, have the same direction as the special solutions x,* throughout the time 
interval, differing from them only in magnitude. 

Now, assemble the e<(f) vectors by columns into the matrix £ (t). The matrix analog of (14) is 

£ = A£-£J(t) } (16) 

where J(i) is the diagonal matrix whose entries are the <r t (t). This is a relationship which is very familiar 
from time-periodic systems. 

We wish to use the e,(t) vectors as the coordinate vectors for describing the solution to the variational 
equations. To this end, define new coordinates y on the tangent space as 

x(t) = £(t)y(t), (17) 

Since £(t) is a nonsingular matrix function of time, at least for to < < < t / , all stability information resides 
within the y variables. Again differentiating (17) and substituting into the variational equations (3) we have 

y = U~ x A£-e-H}y. (18) 

But using (16), this easily reduces to 

y = J(t)y. (19) 

So, this transformation takes the variational equations (3), and replaces them with a set of decoupled , time- 
dependent coefficient differential equations for the variables y, and another set of linear equations (16) for 
the coordinate vectors e(f). We will refer to y as the modal variables for the system, and £(t) as the modal 
matrix. 

The transformation (17) will be canonical if £ is a symplectic matrix for all time. It is possible to so 
normalize £(t 0 ) at the initial time, Siegel and Moser [5], Wiesel and Pohlen [7]. Then £ will stay symplectic 
if £ t Z£ = Z for all time. Taking a time derivative of this and substituting from (16) gives 

£ t A t Z£ + £ t ZA£ - J T ( £ t Z£ ) - ( £ T Z£ ) J = 0. (20) 

Assuming that the modal matrix is at the moment symplectic replaces the quantities in parentheses with Z. 
For Hamiltonian systems the matrix A = Zd 2 H/dX 2 , and since d 2 H/8X 2 is symmetric, the above reduces 
to 

- J T Z- ZJ = 0. (21) 

Simple calculation will show that this is an identity for any diagonal matrix J, so the modal transformation 
(17) is a canonical transformation if £ (<o) is symplectically normalized. The modal equations of motion (19) 
then come from a modal variational Hamiltonian 

K(y) = \y T Z T Jy. (22) 

This is, of course, only a local approximation, ignoring cubic and higher order terms. 
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3 The Voyager 1 Trajectory 

As an example, we have chosen to study an approximation to the Voyager I trajectory from earth past 
Jupiter, and onward almost to Saturn. The trajectory is only approximate, since it was constructed in the 
restricted problem of three bodies. But our method is general, and easily extends to more complex and 
realistic dynamics. The Hamiltonian function is 

h = \(pI + pI+p 2 .) -p*v+p v *-^^-- ^ ( 23 ) 


where 


r 3 = (x - /i) 3 + y 3 + z 1 , 

r\ = (x + l-/i) 3 + y 3 +z 3 . (24) 

Using Voyager I’s known distance of closest approach to Jupiter of 780,000 km and the flight time of 544 
days from the earth to Jupiter, a boundary value problem was posed: starting at Jupiter the trajectory was 
propagated backwards, and at “launch”, it should be 1 A.U. from the sun, and moving tangentially. The 
initial conditions for this trajectory are listed in Table I. In addition, we give a point about 100 days prior 
to close approach. 


Table I 



Launch Initial Conditions 


0* 

Pi 

x y 

+3.762779457438691 x 10 _J +1.886728183030001 x 10” 1 

-2.991922520851858 x 10° +5.825188979188912 x 10" 1 

z 

0.0 

0.0 


Intermediate Initial Conditions 



x y 

z 

0f 

-8.775683982224044 x 10 -1 -4.272485353294678 x 10 -2 

0.0 

Pi 

-8.227825491293955 x 10 -1 -7.044554120116425 x IQ -1 

0.0 


Over this trajectory, a final time of 1.5 dimensionless time units will take the spacecraft from the earth 
almost to the orbit of Saturn. This is shown in Fig. 1, in the rotating reference frame usually used for the 
restricted problem. The regional Lyapunov exponents are ±5.637879, ±5.08635, and ±2.251739. Examining 
the eigenvectors, the first and third are modes in the orbital plane, while the second is purely an out of plane 
mode. This leads to amplification / contraction of initial errors by multiplicative factors of 4706, 2058, and 
29.3. Over an infinite time interval we would expect two zero Lyapunov exponents, and the third mode is 
much the least dramatically unstable of the three. It probably corresponds to an initial displacement along 
the trajectory itself. 

To gain further insight, adjacent trajectories have been examined in the modal space. Integrating the 
nominal trajectory and a nearby orbit, the difference x = X(t) - X 0 (<) was converted into modal variables 
with y = £~ l x. Initial conditions were chosen to excite only one mode at a time, and to explore the limit of 
the linearization inherent in our solution. A shorter arc of the trajectory, spanning ±100 days before and after 
flyby was chosen in order to avoid the extreme differences between initial and final modal amplitudes. Initial 
conditions for this arc are also given in Table I, and the limits of this portion of the trajectory are indicated 
in Fig. 1. Over this interval, the first two modes expand/decay by a factor of about 400, while the third mode 
is nearly static with an expansion/contraction factor of 1.32. The corresponding Lyapunov exponents are 
±20.213590, ±20.161805, ±0.935968. The study of a shorter, less violently expanding / contracting interval 
makes it possible to see the entire modal behavior on graphs. It also emphasizes that neither the Lyapunov 
exponents nor the modal vectors £ are invariant to changes in the trajectory arc studied. 

Fig. 2 shows the behavior of the Lyapunov exponents through this time interval. That is, the figure 
plots the running values of A,(t) starting 100 days before Jupiter approach, and ending 100 days afterwards. 
The final values are the Lyapunov exponents used to decouple the entire trajectory arc, but intermediate 
values show where error growth occurs. Obviously, the immediate vicinity of the close approach is a time of 
explosive error growth. But after the flyby some of this error growth decays again. 
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Figure 3: Modal trajectories on the tangent space for mode 1, with Lyapunov exponents ±20.21359. 

Fig. 3 shows the behavior of the first mode over the ±100 day arc. The modal vectors for this mode lie 
entirely within the orbital plane of Jupiter. Since the modal equations of motion (19) are time dependent, 
Fig. 3 is not a true phase portrait. Rather, initial conditions were scaled by a constant factor to find where 
the system visibly departs from linearity. Since y x grows by a factor of about 400 in this interval, while 
V< shrinks by the same factor, initial conditions are virtually on the vertical y 4 axis, while all trajectories 
terminate nearly on the horizontal y x axis. The linear regime appears at the core of the figure as a symmetric 
region where trajectories scale linearly. Most of the loop is traversed in a very short period of time about 
closest approach. Over most of the time interval trajectories are slowly departing from the vicinity of the y 4 
axis before flyby, or converging towards the y x axis after closest approach. These trajectories were calculated 
with initial values of ya, ya, ys, and y$ zero. These values stayed zero, confirming the success of the modal 
transformation. 

Fig. 4 shows tangent space trajectories for the second mode. This mode lies entirely along the z, p 9 
directions in phase space, and the error growth / shrinkage factor is again about 400 over this trajectory 
arc. The truly linear regime again appears at the core of the figure, while the outer trajectories show visible 
departures from linearity. Also like mode 1, virtually all of the outer loops are traversed in a short time 
interval bracketing closest approach. 

Mode 3 appears to span the in-track direction and the normal vector to the constant Hamiltonian surface. 
Since errors in these two directions are, in the long run, nearly static, the author suspects that if extended to 
infinity that this mode would generate the predicted pair of zero Lyapunov exponents. (A pair since in any 
autonomous system one Lyapunov exponent must be zero, and as a canonical system its Lyapunov exponents 
must occur as positive / negative pairs.) Fig. 5 shows some trajectories for this mode. The expansion / 
contraction of amplitude near close approach is so dramatic for this mode that all the initial and final points 
for this mode are at the origin on this scale plot. Actual initial modal amplitudes were of the order of IQ" 7 , 
and expand briefly by over four orders of magnitude near close approach. The individual trajectories seem 
to leave the origin, and then virtually retrace their outward path in returning. 
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4 Maneuvers 

Since three Lyapunov exponents are positive, while three are negative, we have examined a maneuver strategy 
which attempts to zero the unstable modal amplitudes Vi, Vi, 1/3- The stable modal amplitudes are far less 
important, since they are decreasing any ways. However, we can only perform impulsive changes m velocity in 
the physical space, and we have just stated our maneuver goals in the modal space. Just before a maneuver, 

we have 

x(t~) = S{t)y(t~), ( 25 ) 

while just afterwards we have 

x(r) + *x = £(*)y(r) + £*y. (26) 

Subtracting the two equations above produces 6x = €6y. Now, physically we must have Sx T = (° T > 
since an impulsive maneuver cannot change the position vector, and since the momenta are re y t e 
components of the inertial velocity on the rotating frame axes. In partitioned form, then, the maneuver 

conditions become f x f c 

En E\2 


f 0 \ _( Su e n \ ( \ 

V Sv ) \ £i\ £n J v 6 y*-« ) ’ 

in three by three vector partitions. To zero the unstable modal amplitudes we must have 

SyTs = (-i/i, -i/a, -vs)- 


(27) 


(28) 


The changes in the stable modal amplitudes are not within our control. The first three rows can be solved 
to yield 


<5y 4 - 


( -Vl \ 

i = SxtSn I -V2 I • 


Then, the second three rows give the desired maneuver as 

Sv = {£21 + } I — 



(29) 


(30) 


An alternate strategy is suggested by the fact that only two of the modes experience significant expansion, 
while mode three is nearly static. Attempting to zero only the amplitudes of modes one and two enables us 
to minimize the velocity change required in the maneuver. Rewrite (27) as 


( 5yi-3 ^ _ f f i_2* \ ( 0 V 

\ Sy+-6 ) \ £21 ^22 J \ / 


(31) 


(The £-■ 1 are three by three blocks of the inverse matrix £~ l .) Then, we wish to minimize At^ subject 
to the constraints nAv = -yi and f 3 Av = -y 2 , where the c,- are the first and second rows of £jj . Using 
two Lagrange multipliers \ it the minimum amplitude maneuver is given by the solution to the five linear 

equations 


2Av{ + Ai«ii + — 0, i — 1 , 2 , 3 , 

eiAv = —in, 
e 2 Av = — 1/2. 


(32) 


To study this maneuver strategy, we have begun at 100 days before flyby with a unit error in either yi 
or 1 / 2 - Of course, this is a linear problem, and scales linearly to other initial errors. Then the error was 
corrected at maneuver time t m , and the new state, including nonzero amplitudes m the modal variables y< 6 , 
was propagated to 100 days after the flyby. The results for yx are shown m Figs. 6 and 7. Fig. 6 shows the 
required Av maneuver amplitude as a function of the maneuver time m order to eliminate a unit error m yi 
at the start of the trajectory arc. As expected, the error is much cheaper to correct early in the trajectory, 
and becomes very expensive to correct after flyby. Fig. 7 shows the final values y 4 - s(</) as a function of the 
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&rr^" itUd ' *» »'«« » Utnl unit t 

tm ' The error “ vastl y cheaper to correct before flyby. 


the final time) if the initial yi error is corrected my tii^eV^? 1 ^ Unwanted “crease (again, al 

Vu y *' remaina “sensitive to Vl corrections before close a^mad,* 1 ^ aPPr0ach 1116 col U u S a te variable tc 
early mmenver has considerable time in which to decay This t n tT* mtroduced “to V* by an 

V*(t } ) shows a Unear growth with t m for maneuvers nerfoTm^H ft 5“ ^ the flyby - 40(1 the finaI e »or 
of plane mode y 6 is decoupled from the planar mode* «n * d th * tune of clo8est a PProach. The out 
Figs. 8 and 9 show the analogous resuSToJ ° T*** * * “ itial value <* ’«o. 

cost of correcting an out of plane error soars enormowly jusTat 71“ “ 8tabIe vertical mode The 
drops to a lower, nearly static value. The nearly stitfc ** T P f ?° 8est a PP™ach, and afterwards 

mode growth / contraction in this problem occms ^^ds » due to the fact that most of the 

amplitudes themselves become almost static. Fig 9 shows the* S m Af f. er ^ ards moat of the modal 
maneuver tune. Since the z mode is decoupled from boTh „! a “P^tudes as a function of the 
mitial zero values. However, the stable mode y, comugle toX n0t CXcited from their 

maneuvers performed after the close approach As 1,1! h unstable vertical mode y, is excited by 

significantly decay if the maneuver is performed before the ttthU'h th “ “ due to the fact that W can 
is performed after flyby. P Cf ° re the flyby - but does not greatly decay if the maneuver 


5 Discussion and Conclusions 


cttpTo^^ separation for thne dependent Unear systen* 

—• b , m .p Pr L mat ”n: t sail e** ; h ? ~ m beh *™' - 

make it possible to assess the effects of initial errors i Ih77 of Ju P lter - The modal separation variables 
to the actual size of the errors, and making maximum ?!? take “ J° correct them without reference 
modal transformation offers. tbe P°® 8, ble dynamical decoupling that the 
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Figure 9: Final modal amplitudes yi{tj) as a function of the maneuver time t m used to correct an initial 
unit error in the unstable vertical mode 


a canonical transformation might be used to construct a perturbation theory about the underlying reference 
orbit, including higher order terms in the modal Hamiltonian (22) as the perturbation source. This might 
significantly extend the region of validity of the linearization of the trajectory. 
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